Simulating Physical System with MATLAB
This is a preliminary guide to simulate a dynamical system using MATLAB. At the end, the audience should be able to simulate 2-DOF manipulator with PD control. But first we will look at the first example, pendulum. Keep in mind that there might be other more efficient ways to simulate the same dynamical system. However, this serves as a simple example for those who are interested in making a quick simulation of their systems. There are several things to keep in mind for writing the code for simulation. # How long does the simulation is going to take ? # What is the smallest time step of the simulation ? Are they fixed ? # Do you have the state-space representation of your system ? # What is the initial condition of your system ? # Is your system a closed-loop system ? any inputs ? # How will you interpret your result ? animation, or graph ? MATLAB Environment and Workspace One can write the entire simulation in one single script as well as multiple functions. For the sake of simplicity and debugging purposes, we will use one single script to simulate the system. The format of the code that I use is the following. # Setting up parameters and initial condition # Defining the dynamics of the system # Using ODE solver to find solutions # Interpret the solution with graph and animation These are criteria for the simulation of the pendulum. * The initial time is zero. The duration of simulation is 10 seconds. * The initial condition is when the pendulum perpendicular to the direction of gravity. The initial motion is at rest. * The pendulum is assumed to be point mass. The length of the pendulum is 1 meter. And the mass is 1 kilogram. The gravitational acceleration is 9.81 meter per second square. And the damping coefficient is 0.5 kilogram per second. * The dynamics of the pendulum can be described with the following state-space representation. \left[ \begin{array}{c} \dot{x_1} \\ \dot{x_2} \end{array} \right] = \left[ \begin{array}{c} x_2 \\ -\frac{g}{l}\sin(x_1)-\frac{b}{ml^2}x_2 \end{array} \right] where x_1 is the angular displacement in rad x_2 is the angular velocity in \frac{rad}{s} m is the mass of the pendulum in kg l is the length of the pendulum in m g is the gravitational acceleration in \frac{m}{s^2} b is the damping coefficient in \frac{kg}{s} Setting Up Parameters and Initial Condition In case of testing the simulation multiple times with different value of parameters, it is more convenient to declare variables in the beginning of the script. These variables will store in MATLAB's workspace. One can access, change, and clear these variable at any time in the code as long as they are declared before. For the sake of consistency, we will define parameters into 3 different categories. # Time Elements # Physical Parameters # Initial Condition Time Elements One needs to know how long is the simulation going to last and when the simulation will start. Usually, the start time of most simulation is default as 0. In the code, we will use "startTime" for start time. Instead of end time, we will use duration to describe how long that the simulation will last. We will also denote the duration of the simulation by "duration"' '''Because ODE solver is a numerical tool to solve differential equations, the solution can be computed with discrete time step. This '''time step' can be either fixed or variable. If one decide to use fixed time step, a sample rate for the time step has to be defined. The sample rate in this context refers to how many sample does it has in 1-second period. We will denote sample rate and fixed time step by "fs" and "dt", respectively. Finally, we need to generate a time vector based on the given information if one decides to use fixed time-step. Time vector is a numerical array that contains time-stamp with the specified interval. We will denote the time vector by "t". If one decides to use variable time step, generating time vector can be ignore in this step. The following is the sample code in MATLAB. startTime = 0; duration = 10; fs = 200; % Hz dt = 1/fs; t = startTime:dt:startTime+duration; If one knows the time interval of the simulation, the following code also works. startTime = 0; endTime = 10; fs = 200; % Hz dt = 1/fs; t = startTime:dt:endTime; Physical Parameters These quantities might occur at different places in the code. Therefore, it is convenient to declare variables that represent all parameters in the workspace. The following is the sample code in MATLAB for declaring physical parameters of the system. m = 1; l = 1; g = 9.81; b = 0.5; Initial Conditions It is essential to know the initial states of the system. It is the same idea of analytically solving ordinary differential equation. One needs to know the initial conditions of the differential equation. The total of n initial conditions are required for n^{th} -order system. We will also store the initial conditions in a row vector. The order of the element in the vector is important. The order has to be the same with the state space. In this case, the first initial condition is the initial angular displacement of \frac{\pi}{2} . The last initial condition is zero because the motion starts at rest. The following is the sample code in MATLAB for assigning the initial conditions of the system. We will denote the initial condition by "IC". IC = [ pi/2 , 0 ]; Defining the dynamics of the system Since we are using only one script, we will introduce the anonymous function using the function handle. Essentially, we will describe the dynamics of the system using a function that can be called by a handle without writing additional MATLAB files. Because we are using ODE solver in MATLAB, the change in state vector have to be assigned as a column vector. Inside the function handle, we will write the rate of change of each state variables. The following is the sample code in MATLAB for describing the dynamics of pendulum using anonymous function. dynamics = @(t,x)(; -(g/l)*sin(x(1))-b*x(2)/(m*l^2)) ; Notice: x refers to a state vector. Therefore, in order to access x_1 , we write x(1) to access the first element in the vector x. Using ODE solver to find Solutions There are varieties of ODE solver in MATLAB. Each solver can handle different types of differential equations. In our case, we will use ode45 to solve our system. As an overview of ode45, this function will take a function handle that describes the differential equations as well as time vector and initial conditions. The function will compute the state variables that corresponds to each time stamp as well as the time vector. The format of ode45 is look like the followings. T,X = ode45(handle,timeVector,initialCondition); If "timeVector" is given as an array of time stamp, the output time vector will be the same. But if the given timeVector is in the form of array of two elements that contains start time and end time, ode45 will generate a variable-time-step vector. Let say that ode45 generates the total of 500 samples. The data structure of T is double (500x1). X is the value of state variables at every given time stamp. The data structure of X is double (500x2) for 2 state variables. In order to access the velocity vector, one can use v = x(:,2);. The following is the sample code in MATLAB for generating solution of state variables using ode45. T,X = ode45(dynamics,t,IC); % one can replace T with ~ Interpret the Solutions with Graph and Animation At this point, you are much done with solving the solution. The rest of the work is to display your result. Therefore, it is important to understand the solutions from ODE solver. One can store state variables in other variables with different names to prevent confusion. We will use "x_1" and "x_2" for the angular displacement and angular velocity at all given time stamp , respectively. Trajectory Graph (against time) One can plot the both result against time to see how the pendulum behaves as the time progresses. We will use the function "plot" for displaying a graph. The following is the sample code in MATLAB for generating a plot with solutions from ode45. x_1 = X(:,1); x_2 = X(:,2); figure(1) % start a new figure #1 hold on; % hold for multiple plots plot(t,x_1,'b'); plot(t,x_2,'r'); legend('x_1','x_2'); xlabel('t s'); ylabel('x_1,x_2'); title('x vs. t'); The entire script can be written as follows. startTime = 0; duration = 10; fs = 200; % Hz dt = 1/fs; t = startTime:dt:startTime+duration; m = 1; l = 1; g = 9.81; b = 0.5; IC = [ pi/2 , 0 ]; dynamics = @(~,x)(; -(g/l)*sin(x(1))-b*x(2)/(m*l^2)) ; ~,X = ode45(dynamics,t,IC); x_1 = X(:,1); x_2 = X(:,2); figure(1) hold on; plot(t,x_1,'b'); plot(t,x_2,'r'); legend('x_1','x_2'); xlabel('t s'); ylabel('x_1,x_2'); title('x vs. t'); Animation (in Cartesian Space) In order to create an animation, it is important to know the general look of the system. Since we are going to create an animation of a pendulum, the plot will rely on the kinematics of the pendulum. In other words, we have to know the geometry of the pendulum. Schematic of Pendulum Based on the geometry of the pendulum, kinematics at the end of the pendulum can be expressed as the following. p_x = l\sin(x_1) p_y = -l\cos(x_1) The idea behind the animation is to display one frame after the other in sequence. In MATLAB, we will set the first frame, which includes a plot that represents the pendulum. After that, we will use for-loop to iterate through each data point and change the setting of the original plot. We will not display every single frame since MATLAB might not be able to handle the speed. We will plot after every 25 frames. A function "pause" is required in order to make the animation run according to the time vector. The following is a sample code for creating animation of the pendulum. hinge = 0; link = -l*cos(x_1(1)); p = plot(link(1),link(2)); axis(l -l l); axis equal; lim(1,:) = xlim; lim(1,:) = (lim(1,2)-lim(1,1))/2; temp = ylim; lim(2,:) = ylim; str = sprintf('Time: $0.2f sec',t(1)); pt = title(str); step = 25; for i = 1+step:step:numel(t), link = -l*cos(x_1(i)); set(p, 'XData', link(1) , 'YData', link(2)); set(pt, 'String', sprintf('Time: %0.2f sec', t(i))); axis(lim(2,:)); pause(t(i)-t(i-step)); end If one would like to display both result, the following is the example code. %% Set up Time Element startTime = 0; duration = 10; fs = 200; % Hz dt = 1/fs; t = startTime:dt:startTime+duration; %% Set up Physical Parameters m = 1; l = 1; g = 9.81; b = 0.5; %% Set Up Initial Conditions IC = [ pi/2 , 0 ]; %% Set Up Dynamics (May require additional file for more complex function) dynamics = @(~,x)(; -(g/l)*sin(x(1))-b*x(2)/(m*l^2)) ; %% Run the solver ~,X = ode45(dynamics,t,IC); %% Display State Trajectory vs. Time x_1 = X(:,1); x_2 = X(:,2); figure('Color','white'); subplot(2,1,1); hold on; plot(t,x_1,'b'); plot(t,x_2,'r'); legend('x_1','x_2'); xlabel('t s'); ylabel('x_1,x_2'); title('x vs. t'); %% Animation subplot(2,1,2); hold on; % Set up first frame hinge = 0; link = -l*cos(x_1(1)); p = plot(link(1),link(2)); axis(l -l l); axis equal; lim(1,:) = xlim; lim(1,:) = (lim(1,2)-lim(1,1))/2; temp = ylim; lim(2,:) = ylim; str = sprintf('Time: $0.2f sec',t(1)); pt = title(str); step = 25; for i = 1+step:step:numel(t), link = -l*cos(x_1(i)); set(p, 'XData', link(1) , 'YData', link(2)); set(pt, 'String', sprintf('Time: %0.2f sec', t(i))); axis(lim(2,:)); pause(t(i)-t(i-step)); end